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Abstract 

The progression of liver fibrosis in response to chronic injury varies considerably among individual patients. The underlying 
genetics is highly complex due to large numbers of potential genes, environmental factors and cell types involved. Here, we 
provide the first toxicogenomic analysis of liver fibrosis induced by carbon tetrachloride in the murine 'genetic reference 
panel' of recombinant inbred BXD lines. Our aim was to define the core of risk genes and gene interaction networks that 
control fibrosis progression. Liver fibrosis phenotypes and gene expression profiles were determined in 35 BXD lines. 
Quantitative trait locus (QTL) analysis identified seven genomic loci influencing fibrosis phenotypes (pQTLs) with genome- 
wide significance on chromosomes 4, 5, 7, 12, and 17. Stepwise refinement was based on expression QTL mapping with 
stringent selection criteria, reducing the number of 1,351 candidate genes located in the pQTLs to a final list of 11 cis- 
regulated genes. Our findings demonstrate that the BXD reference population represents a powerful experimental resource 
for shortlisting the genes within a regulatory network that determine the liver's vulnerability to chronic injury. 
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Introduction 

Liver fibrosis is a common consequence of chronic injury. 
Inducing agents vary from hepatotoxins, metabolic disorders and 
autoimmune reactions to viral infections. A characteristic feature 
of the fibrotic response is the ongoing repair mechanism resulting 
in an excessive accumulation of extracellular matrix [1,2]. Fibrosis 
may progress to liver cirrhosis, which is characterized by severe 
distortion of liver architecture and impaired function. Of note, in 
patients with an exposure to similar environmental risk factors, the 
progression of liver fibrosis varies markedly. Based on the rate of 
fibrosis progression, patients may be classified as 'slow' or 'rapid 
fibrosers' [3]. These interindividual differences have been attrib- 
uted to environmental, but also to genetic (and epigenetic) factors 
[1,4,5]. Several fibrogenic gene variants have been identified, e.g. 
complement component 5 {He) [6] , other chemoattractants such as 
the chemokine CXCL9 [7] and the chemokine receptor CXCR3 [8] 
or metabolic enzymes like the triglyceride hydrolase adiponutrin 
{PNPLA3) [9,10]. Moreover, two recent genome-wide association 
studies [11,12] identified a set of novel potential susceptibility 



genes for liver fibrosis, including PNPLA3, but no specific networks 
underlying fibrogenesis were reported. However, due to the large 
number of factors involved, the systematic identification of genetic 
determinants and networks affecting hepatic fibrosis remains a 
major challenge. 

Systems genetics is a powerful method to dissect the underlying 
mechanisms of complex traits, including predisposing gene 
networks and environmental variants [13,14]. The key experi- 
mental set-up is to make use of a genetic reference population. 
Here, we availed of the BXD set of recombinant inbred (RI) 
mouse lines, which are inbred progeny of F2 intercrosses of the 
inbred mouse strains C57BL/6J and DBA/2J [13,15]. RI lines are 
especially suited as a mapping panel, since they form an 
immortalized set of isogenic lines [16,17], and a large number of 
animals and phenotypes per genome can be analyzed under 
standardized experimental conditions, thus lowering environmen- 
tal noise. This improves the yield of information for the detection 
of genetic loci linked to trait variation, known as quantitative trait 
locus (QTL) mapping [13,15]. In previous studies we have 
demonstrated that the parental strains C57BL/6J and DBA/2J 
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show significant phenotypic variation of key fibrogenic parameters 
and therefore differ in their fibrosis susceptibility [6,7,18]. Since 
these strains also vary in four million genetic sites across their 

genome [19], they provide the phenotypic and genetic diversity 
necessary for mapping studies in liver fibrosis. Furthermore, with 
more than 13,000 genetic markers and over 3,000 phenotypic 
records the BXD lines are one of the best-characterized murine 
reference panels [13,20,21]. 

Our aim was to determine new gene variants that affect hepatic 
fibrosis and to apply a systems genetics approach for the 
identification of gene networks that are critical for fibrosis 
phenotypes. Therefore, we characterized differences in fibrosis 
susceptibility of BXD lines after induction of hver fibrosis with 
carbon tetrachloride (CCI4) and generated toxicogenomic, hepatic 
expression profiles by microarray analyses. Afterwards, we 
associated the genetic variation in our population with transcript 
variation in order to identify determinants of gene expression in 
liver fibrosis. In addition to single QTL and gene-gene interaction 
studies, the combination with expression genetics provided novel 
insights into potential networks modifying hepatic fibrogenesis. 

Methods 

Animals and experimental design 

C57BL/6J, DBA/2J, B6D2 Fl hybrids and BXD lines were 

obtained from The Jackson Laboratory (Bar Harbor, ME) or from 
Oak Ridge Laboratory (lines BXD43, BXD51, BXD61, BXD62, 
BXD65, BXD68, BXD69, BXD73, BXD75, BXD87, BXD90), 
and were bred in the facility of the Neurobsik consortium from the 
VU University Amsterdam. The mice were maintained in a mouse 
facility under controlled environmental conditions. 

In addition to the parental strains and the Fl hybrids, we 
studied 35 BXD lines with an average of six mice per sex and line, 
resulting in a total of 581 mice. Liver fibrosis was induced at eight 
weeks of age. For a period of six weeks, CCI4 was administered by 
intraperitoneal injections twice weekly (0.7 mg CCLi/kg body 
weight in mineral oil, final volume 50 |J.l). Forty-eight hours after 
the last injection, the animals were anaesthetized with isofluran 
and killed by cervical dislocation. Instantly, blood was collected 
from the vena cava inferior, and tissue samples of liver and spleen 
were harvested. Liver samples were divided into five separate 
lobes. Whole liver weight, spleen and body weights were noted. 

The animal studies were conducted according to all relevant 
welfare regulations and the Animal Care and Use Committee for 
Saarland University approved the protocols (TV Nr. 10/2008). 

Phenotypic cliaracterization of iiepatic fibrosis 

We measured the following quantitative CCli-induced pheno- 
types: hepatic collagen contents (hydroxyproline levels and 
collagen areas) as quantitative measures and fibrosis stage as 
semiquantitative measure in histological liver sections. In 35 BXD 
lines, C57BL/6J, DBA/2J and Fl hybrids, collagen contents were 
determined in liver hydrolysates from snap frozen specimens of the 
right hepatic lobe. The assay is based on photometric measure- 
ment of the collagen specific amino acid hydroxyproline (Hyp) and 
follows the slighdy revised protocol of JamaU et al. [18,22]. 

For the histological assessment of liver injury, formalin-fixed left 
lobes (4%, v/v) were available from 29 BXD lines, strains C57BL/ 
6J and DBA/2J as well as B6D2 Fl hybrids. Each lobe was cut 
into 3—4: cross sections and embedded in one paraffin block. To 
detect collagen fibers, paraffin sections were stained with Sirius red 
[18]. The staging of fibrosis was performed using a semi- 
quantitative scoring system adapted from the system of Batts 



and Ludwig [18,23,24], principally difierentiating the stages FO to 
F4 (T-scores'). 

Furthermore, stained collagen areas were quantified by 
morphometric analysis, using a semiautomatic system for image 
analysis (Stingray F146C IRF Medical camera, Vi" type progres- 
sive scan CCD, Germany, and HistoQuant image morphometry 
software, SDHistech, Budapest, Ungary). Mean collagenous areas 
((xm^) were calculated by setting a threshold capturing Sirius red 
stained areas of collagen. One representative field (magnification 
lOOx) was chosen from each liver section (avoiding arteries of a 
diameter >100 [tm), and the mean percentage of the stained area 
to whole area (field of vision) was calculated. Liver injury was 
assessed by serum alanine aminotransferase (ALT) activity. After 
CCI4 challenge, blood was collected in a terminal procedure as 
described above. Blood was centrifuged for 20 min with 2000 at 
4°C. Serum was diluted with 0.9% (v/v) NaCl, and ALT levels 
were determined in the central laboratory of Saarland University 
Medical Center according to the IFCC reference method (Cobas, 
Roche Hitachi, Indianapohs, IN) [25]. 

Microarray analysis of hepatic expression profiles 

Total RNA was isolated from snap frozen individual liver 
samples (~30 mg) of 30 BXD strains, the parental strains and 
B6D2 Fl hybrids, using the RNeasy mini kit (Qiagen, HUden, 
Germany). Three female mice per strain were analyzed after CCI4 
treatment for six weeks as described above, resulting in a total of 
99 liver samples. 

RNA quality was verified by measurement of the RNA integrity 
number (2100 Bioanalyzer, Agilent, Santa Clara, CA). Whole 
genome profiles of the fibrotic livers were performed using Gene 
Chip Mouse Gene 1.0 ST arrays (Affymetrix, Santa Clara, CA). 
For the normalization of robust multi-array average (RMA) 
intensity estimates of each transcript, RMAs were transformed into 
log2-values. Then the data of each single array was converted to 
their Z-scores, so that each array has a mean of 0 and a standard 
deviation of 1. Subsequentiy all values were multiplied by 2, and a 
value of 8 was added. Accordingly, all final values are positive and 
all datasets have an average of 8 units. One unit of expression on 
this scale resembles approximately a two-fold difference in 
expression. 

The Affymetrix expression data set comprised a total of 34,760 
records, which were assigned to 22,349 annotated genes. We 
determined the mean gene expression values for each probe set. As 
internal control, the strain distribution patterns of eQTLs with a 
Mendelian (monogenic) expression pattern (e.g. Alad, delta- 
aminolevulinate dehydratase; He, hemolytic complement; Tceanc2, 
transcription elongation factor A) were determined to show a 
perfect match to those of their closest markers, verifying that there 
were no errors of strain assignment in this data set (not shown). 

AU normalized transcript data are available in the GeneNet- 
work database (accession ID GN325, database name SUH BXD 
Li\-c-r CC14-treated Affy Mouse Gene 1.0 ST (Junll) R]\L\). 
GeneNetwork is an open-access database that collates genomic 
information of diverse experimental crosses and reference panels 
as well as phenotypic data from miscellaneous research groups 
[26]. 

Statistics 

Data generation, statistical analysis and graph creation were 
performed with SPSS Statistics 2 1 (IBM, Ehningen, Germany). As 
appropriate, mean and median values were further used for QTL 
analysis. Phenotypic robustness for each strain was assessed by the 
standard errors of the means. Mean and median values were 
trimmed by identifying and omitting outUers after graphical 
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inspection of the data in box plots. Trait values of each BXD line 
were analyzed for each sex separately, as well as for the combined 
data sets of female and male mice. All fibrosis trait data were 
uploaded into the GeneNetwork database (accession IDs 14355- 
14396). Pearson's correlation was used to correlate fibrosis data 
among themselves and to BXD phenotypes. 

HeritabUity (h ) of the fibrosis traits was calculated as the ratios 
of between-strain variances to within-strain variances in sex- 
specific data sets [27]. Within-strain and between-strain variances 
were calculated with analy.sis of variance (ANOVA) [28]. 

pQTL analysis. The identification and mapping of pheno- 
typic QTLs (pQTLs) was performed by linking trait data to 
genotypes at known genetic marker loci [13,29]. All phenotypic 
data were integrated into the GeneNetwork database. For the 
identification of single QTLs, interval mapping analyses were 
performed across all chromosomes [30]. The parental strains were 
included in all mapping analysis. For composite interval mapping 
(CIM), a single genetic marker representing an identified QTL 
region was included as covariate, increasing the power to identify 
QTLs on other chromosomes by removing the effect of the pre- 
eminent QTL [31]. CIM was performed for every phenotype, 
choosing genetic markers with highest likelihood ratio statistics 
(LRS) at each single QTL region or interacting loci from pairwise 
interaction scans. 

The significance of a hypotheticcd QTL was estimated from the 
LRS [30]. Genome-wide significance was evaluated by testing 
2,000 permutations [32], which specified a significant threshold 
corresponding to a genome-wide p-value (po) of 0.05, and a 
suggestive threshold corresponding to po — 0.63. Confidence 
intervals of chromosomal regions spanning QTL positions were 
specified as 1 .5 logarithm of the odds (LOD) support intervals. 

eQTL analysis. In this study, we aimed to infer causal 
mechanisms for phenotypic variation. Observing effects on gene 
expression that result from variants in the identified genomic 
region increases the confidence that this locus harbors causal 
candidates underlying the phenotype. Therefore we followed a 
pre-defined selection strategy for candidate genes, herein we chose 
to restrict our candidate search to the genes located in significant 
and phenot)pe-overlapping pQTL regions (Figure SI in File SI). 
Genes were identified using the GeneNetwork/UCSC Genome 
browser [33]. Since gene expression is, at least in part, heritable, 
differences in mRNA expression levels over the panel of BXD lines 
can be used to map regulatory expression quantitative trait loci 
(eQTLs). Therefore hepatic expression levels were used as 
quantitative traits and implemented into interval mapping analyses 
to identify regulatory loci. An expression quantitative trait locus 
(QTL) for a specific transcript was denoted cisQTL if the 
associated marker was localized within a 10 Mb distance of the 
gene position [34]. The respective gene was called a cis-regulated 
quantitative trait gene (cisQTG). 

Integrated search for candidate genes. The lists of 
cisQTGs for the fibrosis phenotypes were refined following an 
explorative in silico data analysis (Figure S2 in File SI). We applied 
the following selection criteria, of which at least one had to be 
fulfilled by the a'jQTG to qualify as a potential fibrogenic 
susceptibility gene: (i) significant Pearson's correlation coefficient 
(p<0.05 corresponding to r>0.36 or r< — 0.36) between a'rQTG 
and any of the fibrosis phenotypes (collagen area. Hyp, F-score); 
(ii) non-synonymous single nucleotide polymorphism (nsSNP) that 
differs between the parental strains B6 and D2. nsSNPs leading to 
amino acid substitutions in the coding regions were identified 
using the GeneNetwork variant browser; and (iii) differential 
hepatic regulation of cijQTGs in expression data sets after CCI4 
challenge and saline treated control livers. QTL regions were 



matched using QTLminer [35] and the GeneNetwork dataset for 
sahne treated hvers in females (accession ID GN312, database 
name GenEx BXD Sal Liver Affy M430 2.0 (Feb 11) RMA 
Females). 

Gene network construction. Pairwise correlation estimates 
of the expression values of designated candidate genes were 
calculated as Pearson's correlation coefficient and presented in a 
circular network graph using Cytoscape software [36]. For 
purposes of illustration the intra-chromosomal correlations of 
genes on the same chromosome were omitted, and only inter- 
chromosomal correlations are shown in the network. Edge colors 
represent positive or negative correlations with r>0.36 or 
r< — 0.36, and node sizes indicate the connectivity of the genes. 

All candidate genes including the three fibrosis phenotypes were 
additionally illustrated in a QTL heatmap. This heatmap 
visualizes the p-values of regulator}' loci identified liy genome- 
wide linkage analysis of the traits, computed on the basis of 
permuation tests (n= 1,000). 

Data access 

All trait data were uploaded into the GeneNetwork database 
"BXD Published Phenotypes Database" (http://www. 
genenetwork.org/): Fibrosis phenotype data (CCI4 treated livers) 
accession IDs 14355-14396; expression dataset (CCI4 treated 
livers, females) accession ID GN325, database name SUH BXD 
liver Affy Mouse Gene 1.0 ST (Junll) RMA; expression dataset 
(saline treated livers, females) accession ID GN312, database name 
GenEx BXD Sal liver Affy M430 2.0 (Feb 11) RMA Females, 
provided by courtesy of Dr. Robert Rooney, Genome Explora- 
tions (Memphis, TN). 

Results 

Liver fibrosis in parental strains 

A significant (p<0.0001) accumulation of collagen was observed 
in all livers after six weeks of CCI4 challenge in comparison to 
untreated mice. The mortality rate of mice after CCI4 challenge 
was 8%. The parental inbred strain DBA/2J is more susceptible to 
liver fibrosis than the C57BL/6J strain, indicated herein by 
significantiy (p<0.05) increased collagen areas and higher hepatic 
collagen (Hyp) contents; this is paralleled by liver injury as assessed 
by serum ALT activities (Figure 1). 

Liver fibrosis phenotypes in the BXD reference panel 

We noted marked differences in liver fibrosis among the BXD 
lines. Hepatic collagen contents varied widely (mean ± SD 
386. 9± 141.5 ng H)p/g liver; Figure 2A), and histopathological 
fibrosis scores correlated significantly (p<0.001) with hepatic 
collagen levels (Figure 2B). Furthermore, we observed significant 
(p<0.0001) line differences with respect to hepatic collagen levels 
and fibrosis scores (Figures 2C-E) as well as clinical-chemical 
parameters (ALT) of liver damage (Figure 2F). Semiquantitative 
fibrosis scores ranged from Fl (perivenular fibrosis) to F4 
(pronounced fibrosis) (Figure 2G— K). 

Heritability of liver fibrosis as determined by mean h^ values 
was similar for all fibrosis traits, ranging from 0.51±0.18 (Hyp) to 
0.57±0.02 (F-score) and 0.59±0.01 (collagen area). Overall, 
hepatic collag(;n ar(;a was identified as tin; most heritable trait 
(h — 0.59). h for hepatic Hyp concentrations showed a difference 
between male (0.36) and female mice (0.87), whereas h^ for the 
other traits did not differ between sexes. 



PLOS ONE I www.plosone.org 



3 



February 2014 | Volume 9 | Issue 2 | e89279 



Systems Genetics of Hepatic Fibrosis 





E 








DBA/2J 





Figure 1. Phenotypic characterization of parental strains after six weeks of CCI4 injections. Liver fibrosis was assessed by morphometric 
(A) and biochemical (B) measurement of hepatic collagen (Hyp) contents. Hepatic inflammation was measured by serum ALT activities (C). Sirius red 
staining of hepatic collagen showed circumferential fibrosis in C57BL/6J mice (D) and pronounced fibrosis in DBA/2J mice (E), corresponding to mean 
F-scores of 2.0±0.1 and 3.9±0.1, respectively. 
doi:1 0.1 371 /journal.pone.0089279.g001 



Genome-wide mapping of liver fibrosis phenotypes 
(pQTLs) 

Single QTL genome scans identified 28 trait associated loci 
(Table SI in File S2) with LRS scores above the suggestive 
threshold (see Methods] affecting liver fibrosis phenotypes. Among 
these significantiy linked loci were detected by composite interval 
mapping (pG<0.05) on chromosomes 4, 5, 7, 12, and 17 (Table 1, 



Figure S3 in File SI). For all loci except the QTL on chromosome 
12, alleles of the fibrosis-susceptible strain DBA/2J increased the 
trait values (Table SI in File S2). Five QTLs on chromosomes 2, 5, 
7, 13 and 15 conferred susceptibility to more than one phenotype 
(Table S2 in File S2), whereas the loci on chromosomes 4, 12 and 
17 were specific for a single phenotype. 
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Figure 2. Phenotypic characterization of the BXD reference panel after six weeks of CCI4 injections. (A) Histogram illustrating the 
distribution of hepatic collagen contents in BXD mice. The mean hepatic collagen concentrations (± SEjare indicated by vertical lines (solid line: BXD 
recombinant inbred lines, 386.9±5.9 ng Hyp/g liver; dashed line: C57BL/6J (B6) inbred strain, 419.3±22.4 |ig/g; dash-dot line: DBA/2J (D2) inbred 
strain, 570.9±25.9 ng/g). (B) iVlean hepatic collagen contents (± SE) stratified according to F-score categories (F1-F4). (C-E) Strain specific mean (± 
SE) phenotype values compared to the overall means of all mice phenotyped (29-35 BXD lines, strains C57BL/6J and DBA/2J, B6D2 Fl hybrids), which 
are represented by the horizontal lines. (C) Hepatic collagen areas (mean ± SE: 2.4±0.1%); (D) hepatic collagen (Hyp) contents (386.9±5.9 |ig/g); (E) 
fibrosis stages (F-scores, 2.3±0.1); (F) serum ALT activities (419.4±25.3 U/l). (G)-(K) Representative liver sections of diverse BXD strains after CCI4 
treatment, illustrating increasing F-scores from FO to F4. (G) FO: no fibrosis; (H) Fl: perivenular fibrosis, initially forming collagen bridges; (I) F2: 
circumferential venular fibrosis w\th incomplete bridging; (J) F3: centro-central fibrosis with complete bridging; (K) F4: pronounced fibrosis with 
complete and broadened collagen bridges. 
doi:1 0.1 371 /journal.pone.0089279.g002 
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We found that 16 of 28 QTLs (46%) were sex-specific, i.e. they 
were only found in data from male or female mice (Table SI in 
File S2). The remaining 1 1 loci were detected in the combined 
datasets; but showed significant effects either in male or female, 
consistent with a predominant phenotypic eflFect of one sex. In 
addition, the chromosome 2 QTL at 174.5-181.5 Mb was 
detected in all datasets tested. QTLs for hepatic Hyp levels were 
mainly based on female datasets; this was in line with the higher h^ 
scores for this trait in female mice. For collagen area and F-score 
QTLs, no predominance of a single sex was observed. 

Genome-wide mapping of fibrosis-associated eQTLs 

Nine pQTL regions on chromosomes 2, 4, 5, 7, 12, 13, 15 and 
17 were further dissected using eQTL mapping (Figure 3). The 
nine selected loci either conferred significant linkage or were 
associated to more than one phenotype. eQTLs were ci.s- 
regulatory loci (ciiQTLs) or franj-regulatory loci (franjQTLs) (see 
Methods). By mapping m-regulatory eQTLs within the nine most 
significant pQTL regions, we identified fibrosis-associated expres- 
sion patterns that were locally regulated within the pQTLs. Table 1 
summarizes the results of the eQTL analysis: On average, the 
pQTLs spanned an interval of 15.2 Mb, and these pQTL regions 
contained a total number of 1,351 annotated genes. The highest 
LRS score was observed for the QTL on chromosome 12, the only 
locus for which fibrosis susceptibility was conferred by the C57BL/ 
6J allele. The QTL on chromosome 7 was the largest QTL with a 
high gene density. Overall, we identified 68 regulatory markers 
(eQTLs) within the pQTL regions. Thirty regulatory markers were 
identified as cuQTLs. Using the upper hmit of suggestive 
thresholds for genome-wide significance as determined by 
permutation tests (LRS&12.0), these markers were linked to 85 
genes in close proximity (<10 Mb). The associated genes are 
potentially cM-regulated genes (n'jQTGs) in the pQTL regions. Of 
note, all cwQTLs also demonstrate frflK,?-regulation of additional 
genes outside the pQTL regions (not shown). Applying this 
combined analysis of pQTLs and eQTLs, we reduced the number 
of potential fibrosis candidate genes from 1,351 to 85 (Table 2, 
Figure S2 in File SI). 

Selecting fibrogenic candidate genes 

In further selection steps, we inferred key regulatory candidates 
of fibrosis among the 85 cwQTGs. This strategy refined the list of 
85 cwQTGs to 55 potential profibrogenic candidate genes that 



fulfilled at least one selection criterion (Table S3 in File S2). First, 
Pearson's correlation of mQTG expression with any fibrosis 
phenotype (coUagen area. Hyp, F-score) identified 30 significantiy 
(p<0.05) correlated genes (indicated by dark gray boxes). Next, we 
determined covariance of hepatic mRNA expression patterns in 
unchallenged and CCI4 challenged livers of the BXD reference 
lines, using QTLminer analysis. This analysis indicated that 
eQTLs could be distinguished into (A) fibrosis-specific cuQTLs 
that showed differential regulation between the basal state and 
after fibrosis induction, and (B) fibrosis-independent mQTLs, i.e. 
the genes were cis-regulated in both groups. We speculate that 
differential regulation of class A genes in fibrotic livers identifies 
more relevant modifiers of fibrogenesis. In total, 45 cmQTGs were 
differentially regulated (class A), while 40 cisQTGs were cis- 
regulated in both normal and fibrotic livers (class B). Finally, we 
identified 169 genes in the pQTL regions with nsSNPs that 
segregated between the two parental lines C57BL/6J and DBA/ 
2J. Among the 85 cijQTGs, 22 have nsSNPs in coding regions. 

In summary, 55 candidate genes listed in Table S3 in File S2 
fulfilled at least one criterion and were either significantiy 
correlated to a specific fibrosis phenotype, differentially regulated, 
or contained an amino acid substitution. Only 3 1 genes passed at 
least two criteria, and merely eleven genes fulfilled all three 
criteria: Afm (afamin), Fanl (FANCD2 / FANCI-associated nuclease 1), 
Hsdl7bl4 (hydwxysteroid (17-beta) dehydrogenase 14), Napsa (napsin A 
aspartic peptidase), Noma (nodal modulator 1), Nin (ninein), Susdl (sushi 
domain containing 1), and four members of the kallikrein lb family 
(5, 21, 22 and 26). Most candidate genes are located on 
chromosome 7 (n = 21), followed by chromosomes 5 (11= 11), 15 
(n = 9), 4 (n = 7), 12 (n = 4), 2 (n = 2), and 13 (n = 1). The QTLs on 
chromosome 5 (3.1-20.1 Mb) and chromosome 17 contained no 
cmQTG complying with the criteria. 

Generating a fibrosis gene network 

We generated a circular network graph (Figure 4), using the 
pairwise correlation estimates of hepatic expression levels for the 
55 fibrosis candidate genes (Table S3 in File S2). The nodes 
represent the genes, which were arranged according to their 
chromosomal localizations, and edges show significant (p<0.05) 
inter-chromosomal correlations between gene expression levels. 
This unique network contains a total of 1 1 5 inter-chromosomal 
correlations; for reasons of simplicity the graph does not show the 
93 intra-chromosomal correlations. The genes with the highest 



Table 1. Chromosomal regions of pQTLs with significant genome-wide LRS values determined by single QTL scans and CIIVl. 
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85.1-97.9 


0.451 


female/both 


Hyp 


4 


17.4 


r56254381 - rsl 3477745 


55.1-73.9 


58.957 


female/both 


Hyp 


7 


16.3 


rs3703247 - rs8255275 


52.8-56.7 


56.761 


female/both 


Hyp 


12 


25.0 


rs3716547 - rs13481511 


60.5-73.3 


-77.257 


female/both 


F-score 


7 


20.3 


rs3703247 - rs8255275 


48.2-53.7 


0.562 


male/both 


F-score 


17 


22.0 


rsl 3483077 - rsl 3483081 


64.9-71.1 


0.516 


female 



Abbreviations and definitions: pQTL (chr): chromosomal position of quantitative trait locus; LRS (max): likelihood ratio statistic, maximum association between 
genotype and phenotype variation; SNP (max): single nucleotide polymorphism with maximum LRS in QTL region; 1.5 LOD support interval (Mb): chromosomal 
region in Megabases spanning QTL position; Additive allele effect: estimate of a change in the average phenotype by substitution of one parental allele by another 
at a given marker position; (— ) values indicate an increase of phenotype by C57BL/6J allele, (+) values an increase of phenotype by DBA/2J allele; Dataset: dataset in 
which the QTL was identified; Hyp: hydroxyproline; CIM: composite interval mapping. 
doi:l 0.1 371/journal.pone.0089279.t001 
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Figure 3. QTLs for hepatic fibrosis in tKie BXD murine reference population. The heatmaps represent significant interval mapping results on 
the indicated mouse chromosomes, separately for male and female mice as well as the combined data set (without inclusion of covariates); the QTL 
plots below illustrate composite interval mapping results (with 'background' QTLs as covariates, restricted to significant QTLs or overlapping loci for 
different phenotypes). Color coding of the heatmaps is as follows: Grey/black regions indicate the absence of genotype to phenotype linkage. Blue to 
green regions correspond to suggestive and significant linkage, respectively, with C57BL/6J alleles being associated with higher trait values. Red to 
yellow regions correspond to suggestive and significant linkage, respectively, with an association of DBA/2J alleles with higher values. 
doi:1 0.1 371 /journal.pone.0089279.g003 



inter-chromosomal connectivity are Mcee (n=13), Tnc (n=12), 
Septll (n=ll), and Thap6 (n=10). Additional analyses of our 
transcriptomic dataset showed that Tnc highly correlates with other 
fibrosis-associated genes, in particular with hepatic expression 



levels of coUagens [Colla2: r = 0.76; p = 0.002; ColSal: r = 0.75; 
p = 0.001) and tran.sforming growth factor P 1 [Tgfhl: r = 0.74; 
p = 0.002). 
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In addition, we illustrated tiie regulatory eQTLs of the 
candidate genes and the loci for fibrosis phenotypes (pQTL) in a 
QTL heatmap (Figure S4 in File SI). Regulatory gene clusters on 
chromosomes 2, 4, 5, 7 and 12 co-localized with loci of fibrosis 
phenotypes, as indicated by arrows in the QTL heatmap. The 
large eQTL on chromosome 7 showed two differentially regulated 
gene clusters: The genes located on distal chromosome 7 (63.1— 
74.2 Mb) displayed additional regulatory loci on chromosome 5, 
which were absent for the genes between 51.2 and 53.5 Mb. 
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Discussion 

Previously we identified individual loci that confer genetic 
susceptibility to hepatic fibrosis in different experimental crosses of 
inbred mice [6,7,18]. Here we report the first systems genetics 
analysis of fibrosis in the BXD murine reference panel that allows 
the integration of multiple traits [13]. The genetically mosaic BXD 
inbred lines display significant variation of quantitative fibrosis 
phenotypes, consistent with polygenic inheritance of liver fibrosis 
[1 ,4,5] . By correlating phenotypes and known BXD genotypes in a 
genome-wide QTL analysis, we identified multiple pQTLs, nine 
with genome-wide significance and several with sex-dependent 
effects. Sex-specific differences are observed in various liver 
diseases such as (non-)alcohohc fatty Uver diseases and hemochro- 
matosis and might be due to sex hormone-regulated mechanisms 
or sex-specific gene variants. In addition to the phenotypic 
characterization of fibrosis, we generated a comprehensive 
expression dataset that represents the first genome-wide tran- 
scriptome analysis of hepatic fibrosis in a murine reference panel. 

By stepwise bioinformatic analyses [37], we were able to reduce 
the number of 1,351 genes located in the nine pQTL regions to a 
set of 55 creedal fibrogenic candidate genes (Figure S2 in File SI). 
For this analysis, we focused our search on the genes that are cis- 
regulated during fibrogenesis. To minimize the false discovery 
rate, the genes underwent a subsequent careful explorative 
analysis. They were considered to be relevant for fibrosis when 
their expression levels correlated with the fibrosis traits or showed 
difierential regulation in healthy and fibrotic livers. In addition, we 
screened these ctjQTGs for nsSNPs in the parental strains of the 
BXD panel, which might structurally and/or functionally affect 
protein functions. After all these steps, eleven genes fulfilled all 
selection criteria [Afin, Fan], Hsdl7bl4, KlklbS, mib21, Klklb22, 
Klklb26, Napsa, Nomo, Nin, Susdl). For the majority of these genes 
there is no established connection to hepatic fibrosis or little 
information about their function in liver, although kallikreins exert 
known functions in the activation of inflammation, wound healing, 
and liver regeneration [38]. In particular, we consider Afm as 
interesting candidate with a so far unknown role in hepatic 
fibrogenesis. Afm is a member of the albumin gene family that was 
shown to be differentially regulated by hepatocyte nuclear factors 
la and Ifi in mice [39]. Interestingly it functions as carrier of 
vitamin E [40], which has recently been reported to ameliorate 
liver fibrosis in fatty liver disease in mice and humans [4 1 ,42] . A 
study by Kim et al. [43] revealed that afamin acts as a chemokine 
activating the Akt-signalling cascade, at least in osteoblasts. 
Because expression profiling showed Afm to be differentially 
expressed in mouse liver, this observation suggests similar 
regulatory effects across organs. 

A further proof of principle is that several of the significant QTL 
regions include potential candidates that have previously been 
associated with fibrosis progression, in particular the chemokine 
ligand CxcllO, Mlh2 (nuclear receptor subfamily 1, group H, 
member 2, a.k.a. LXR), and Tnc (tenascin C). CxcllO encodes a 
chemokine that promotes hepatic inflammation by leukocyte 
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Figure 4. Fibrosis network grapKi generated by correlating tKie hepatic expression levels of the candidate genes. The graph presents 
the inter-chromosomal correlations of 51 genes, except for the four candidate genes AdamtsU, Gm9860, Klk22 and Ogfr with exclusively intra- 
chromosomal correlations. Node color and shape illustrate the chromosomal localization of the gene. The size of each node indicates the degree of 
connectivity, with larger nodes having higher number of correlated genes. The edges show Pearson's correlation coefficients (r) as follows: solid 
green lines: r>0.5; solid red lines: r<— 0.5; dotted green lines r>0.36; dotted red line r<— 0.36. 
doi:1 0.1 371 /journal.pone.0089279.g004 



recruitment [44] . Additional in vitro experiments in primary mouse 
hepatocytes detected a time-dependent induction of CxcllO 
expression levels after treatment with Tgfbl, supporting its role 
as profibrogenic candidate gene (R.H., R.L. and F.L.; unpublished 
observations). The nuclear receptor LXR might exert antifibrotic 
effects, since it was shown to reduce hepatic stellate cell activation, 
and therefore inhibit the production of profibrogenic cytokines 
[45]. Tnc is an extracellular matrix glycoprotein expressed by 
hepatic stellate cells and myofibroblasts in liver, where it increases 
cytokine expression and ameliorates leukocyte transmigration 
[46,47]. Whereas Tnc expression is absent in naive livers, it is 
strongly induced during enhanced cell turnover as seen in wound 
repair [48,49]. It also represents an endogenous ligand of T/r^, 
which promotes innate immune responses during fibrogenesis 
[50,51]. Further correlation network analysis of our transcriptome 
datasets showed that Tnc is highly interconnected with other 



fibrosis-associated genes (R.H. and F.L., unpublished observa- 
tions); in particular, it is significantly correlated with the hepatic 
expression levels of coUagens {Colla2: r = 0.76; p = 0.002; CoBaT. 
r = 0.75; p = 0.001) and Tgfbl (r = 0.74; p = 0.002). 

To illustrate the complexity of fibrosis susceptibility, we 
combined the genes and their expression correlations in a large 
'fibrosis network', indicating the inter-chromosomal connectivity 
of fibrogenic genes. In addition, the genetic architecture of co- 
regulated gene networks were visualized by eQTL heatmaps, 
which show that the regulatory gene clusters of fibrosis suscepti- 
bility genes co-localize with pQTLs (Figure S4 in File SI). These 
findings reflect the selection concept for the candidate genes, with 
gene clusters representing networks with a causal relationship to 
fibrosis phenotypes. Gene repression studies (e.g. siRNA) might 
help to further dissect the directionality of the gene effects and 
fibrosis gene networks. 
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Recently genome-wide association studies (GWAS) in patients 
with chronic liver diseases have identified profibrogenic gene 
variants [5] . Genes associated with fibrosis progression in patients 
with chronic HCV infection were MERTK (c-mer proto-oncogene 
tyrosine kinase), KNF7 (ring finger protein 7), and TULPl (tubby 
like protein 1) [12]. In contrast, the genes PJ\TLA3 (adiponutrin), 
GCKR (glucokinase regulator) and TRIBl (tribbles homolog 1) 
were associated with fibrosis phenotypes in non-alcoholic fatty 
Uver disease [52], with PMPLA3 demonstrating the most consistent 
effects across studies [53]. Although our expression dataset shows 
that all murine orthologs are expressed in fibrotic mouse liver 
(mean expression scores 7.2-10.4; R.H. and F.L., unpublished 
observations), solely PnplaS was located in a pQTL on chromo- 
some 15. However, it was apparenfly not «>-regulated in our panel 
and did not display nsSNPs in the parental strains. The lack of 
cross-species validation might be due to the distinct induction of 
fibrosis in our model, the specific charcteristics of the patients 
included in the GWAS (viral hepatitis, fatty liver disease), different 
phenotypic efiect sizes across species, or non-polymorphic regions 
within the BXD panel. 

This study illustrates how advances in the methodologies of 
systems genetics with the use of a murine reference panel lead to 
the identification of a potential disease network for liver fibrosis. 
The BXD lines represent an appropriate reference population 
with phenotypic segregation of fibrosis phenotypes due to genetic 
variation. Our findings indicate that it is essential not to focus on 
single fibrogenic QTLs, but on gene clusters as modifiers of fibrosis 
susceptibility. Although several existing confounders are being 
controlled for within this approach, further developments in gene 
mapping and functional validation will contribute to the transla- 
tion of our experimental findings to patients with Hver fibrosis. 

Supporting Information 

File SI Figure SI, Graphical overview of the experimental setup 
for the integrative analysis of pQTLs and eQTLs in the BXD 
murine reference panel. Abbreviations: BXD, recombinant inbred 
lines based on parental strains C57BL/6J and DBA/2J; CCI4, 
carbon tetrachloride; DNA, deoxyribonucleic acid; eQTL, 
expression quantitative trait locus; F-score, fibrosis score; Hyp, 
hydroxyproline; n, number; pQTL, phenotypic quantitative trait 
locus. Figure S2, Study design and strategy for the selection of 
candidate genes. Genome-wide association studies of CCI4 treated 
BXD lines identified phenotype-associated QTLs (pQTLs). The 
1,351 genes located in significant pQTLs were investigated further 
by eQTL analyses (see Methods). This allowed the differentiation of 
local {cis-) or distant (fra«j-)-regulation of gene expression. Cis- 
regulated genes (cisQTGs) underwent the following three selection 
steps to refine the list of candidate genes: I) eisQTGs with 
significant correlation with fibrosis phenotypes; II) fibrosis-specific 
cijQTGs that show difierential regulation between basal state and 
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